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ABSTRACT 

The supermassive black hole candidate at the center of M87 drives an ultra-relativistic 
jet visible on kiloparsec scales, and its large mass and relative proximity allow for event 
horizon scale imaging with very long baseline interferometry at millimetre wavelengths 
(mm-VLBI). Recently, relativistic magneto- hydrodynamic (MHD) simulations of black 
hole accretion flows have proven capable of launching magnetically-dominated jets. We 
construct time-dependent disc/jet models of the innermost portion of the M87 nucleus 
by performing relativistic radiative transfer calculations from one such simulation. We 
identify two types of models, jet-dominated or disc/jet, that can explain the spectral 
properties of M87, and use them to make predictions for current and future mm-VLBI 
observations. The Gaussian source size for the favored sky orientation and inclination 
from observations of the large-scale jet is 33 — 44/ias (p± 4 — 6 Schwarzschild radii) 
on current mm-VLBI telescopes, very similar to existing observations of Sgr A*. The 
black hole shadow, direct evidence for an event horizon, should be visible in future 
measurements using baselines between Hawaii and Mexico. Both models exhibit vari- 
ability at millimetre wavelengths with factor of ~ 2 amplitudes on year timescales. 
For the low inclination of M87, the counter-jet dominates the event horizon scale 
millimetre wavelength emission from the jet-forming region. 

Key words: accretion, accretion discs — black hole physics — radiative transfer — 
relativity — galaxies: individual (M87) — galaxies: active (M87) — galaxies: jets 



1 INTRODUCTION 

Messier 87 (M87) is a giant elliptical galaxy in the Virgo 
cluster, known for its galaxy-scale, ultra-relativistic jet. 
Very long baseline inte r ferometry (VLBI) images at 7mm 
(ljunor. Biretta fc Livid Il999l : iLv. Walker fe Wrobell 12004 
Walk er et alT T2008) show extended jet structure on milli- 
arcsecond (mas) scales, emanating from an unreso lved bright 
core. Recent images between 7 mm and 1.3 cm (|Hada et al.l 
1201 lh suggest that this core is coincidental with the cen- 
tral black hole candidate (hereafter M87 refers to the black 
hole rather than the galaxy). M87 is also one of the two 
largest black holes on the sky (along with the Galactic cen- 
tre black hole candidate, Sgr A*). Its mass is ~ 6.4 x 1O 9 M 
(|Gebhardt fe Thomad [2009 : iGebhardt et al.1 l201lh . -1600 
times larger than Sgr A*. At a distance of 16 Mpc.the an- 
gular size, 50 oc M/D, is about 4/5 that of Sgr A*Q 

Recent VLBI observations at 1.3 mm have detected 



source structure in Sgr A* on eve nt horizon scales 
(|Doeleman et al.ll2008l : iFish et al.l l201lh . allowing a direct 
comparison betw een observations and black hole accre - 
tion theory (e.g., iDexter et al.ll20lA iBroderick et alJl201lh . 
These observations also have the potential to detect the 
black hole shadow (iBaxdeenl Il973l : iFalcke. Melia fc Agoll 
l2000b IDexter et al]l2010h . which would provide the first di- 
rect evidence for an event horizon. 

M87 is just as, if not more, promising a mm-VLBI tar- 
get as Sgr A*. M87 is in the Northern sky, offering longer 
mutual visibility with current telescopes. Its large black hole 
mass gives a proportionally longer dynamical time, so that 
its event horizon spans —1 light-day. This means that, un- 
like in Sgr A*, Earth- aperture synthesis could be used to 
fill in the u-v plane, potentially (with additional telescopes) 
allowing the creation of an image of the source directly, in 
contrast with the model-dependent Fourier domain fitting 
techniques necessary for Sgr A*. Micro- arcsecond (//as) res- 



* E-mail: jdexter@berkeley.edu 

1 The p revious black hole m ass estimate was a factor of ~ 2 throughout, but discuss the effects of the mass on our models in 

smaller (|Marconi et al.|[l997l V We use the new black hole mass ^6.31 
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olution of mm-VLBI also offers the possibility of imaging the 
jet launching region, which would provide the opportunity 
to compare directly with physical models of jet formation 
in the immediate vicinity of the black hole. The most re- 
cent mm-VLBI c ampaign observed M87 as well as Sgr A* 
(jFish et al.l l2011). but the data are as yet unpublished. 

The spectral properties of M87 are well known: it is an 
inverted radio source with a power law tail extending from 
the spectral peak in the millimetre to the optical. Previ- 
ous semi-analytic work has modeled the low-frequency radio 
emission as arising from synchrotron radiation in either an 
advection-dominated accretion flow ([Reynolds et al.l 1996; 
Di Matteo et al 2003) or a "truncated" accretion disc (Yuan 



2006bh as well as other improvements (iMcKinnevi l2006al : 



2000l : lBroderick fe Loebll2009l ). Both types of models can fit 



the spectrum. Synthetic jet images have also been produced 
as predic tions for mm-VLBI from one of these semi-analytic 
models (jBroderick fe Loebll2009h . 

Maps of synchrotron emission from numerical simula- 
tions of jets have also been compared to observations of 
M87 (Zakamska, Begelman & Blandford 2008; Gracia et al. 
2009), but these jets are input by hand rather being formed 
self-consistently from an accretion flow. Conversely, axisym- 
metric general relativistic MHD (GRMHD) simulations of 
the a ccretion flow have been used to fit the spectrum o f the 
core (jMoscibrodzka et al.ll20lH : iHilburn fe Liandl201lh . 

Jet formation has rec ently become accessible to gl obal, 
3D GRMHD simulations (jMcKinnev fe Blandfor3l2009h . In 
their sim ulation initialised with a dipolar magnetic field 
(MBD in lDexter et~ai1l2010L and here), an ultra-relativistic 
jet is produced self-consistently from an accretion disc and 
propagates stably out to 1000 M before it interacts with 
the simulation boundary^ We perform general relativistic 
radiative transfer via ray tracing to create the first time- 
dependent spectral models of M87 from a GRMHD simu- 
lation, and use the models to make predictions for current 
and future mm-VLBI observations. The simulation and the 
relativistic radiative transfer method used to compute ob- 
servables are described in £[2j along with the details of the 
radiative disc/jet models. Fiducial models for the M87 spec- 
trum are identified in 331 and the resulting millimetre images 
and variability are studied in and JS] The implications 
of our results are discussed in §6\ along with the many un- 
certainties in constructing the models. We summarize the 
results in JJl 



2 METHODS 

2.1 Simulation Data 

iMcKinnev fe Blandfordl {2009) carried out 3D global 
GRMHD simulations starting from a hydrostatic torus 
in whose angular momentum was aligned with the black 
hole spin axis for a spin value a/M = 0.92. The ini- 
tial torus pressure maximum was located at r = 12 
M. The GRMHD code used was a 3D version of HARM 
(|Gammie. McKinnev fe Tothll^Ooi iNoble etalll2006h with 
fourth order interpolation and time-stepping (|McKinnevl 



Tchekhovskov, McKinnev & Naravanll2007n . 

This simulation conserved total energy and neglected 
radiative cooling. Energy- conserving simulations convert the 
energy lost to grid-scale numerical magnetic reconnection 
into heat. This is appropriate for studying radiatively in- 
efficient sources such as Sgr A* or M87, where the accre- 
tion flow is hot and geometrically thick The full azimuthal 
domain was included, but only marginally resolved with a 
resolution of 256x128x32 in (r, The effective reso- 

lution was higher than in fully second order schemes due 
to the high order scheme used. The coordinates were reg- 
ular but warped, with the resolution concentrated toward 
the mid-plane at small radius to resolve the disc and to- 
wards the pole at large radius to resolve the jet. The total 
duration was 3500 M, with roughly constant radial profiles 
in accretion rate and angular momentum out to r c± 10 M 
by t rsj 3000 M. The snapshot used for spectral fitting in 
^3] is from t = 3000 M, after an approximate quasi- steady 
state was established in the inner disc. The final 2000 M 
of the simulation is considered when we study time- variable 
properties of the radiative models. 



2.2 Ray Tracing 

We performed relativistic radiative transfer on the simula- 
tion data via ray tracing using the code grtrans ([Dexter! 
2011). Starting from an observer's camera, rays are traced 
backwards in time toward the black hole assuming they are 
null geodesies (geometric optics approximation) . using the 
public code geokerr described in lDexter &; Agoll (|2009l ). In 
the region where rays intersect the accretion flow, the radia- 
tive t ransfer equation is sol ved along the g eodes ic ([Broderickl 
2006) in the form given in lFuerst fc Wul (|2004h . which then 
represents a pixel of the image. This procedure is repeated 
for many rays to produce an image, and at many time 
steps of the simulation to produce time-dependent images 
(movies) . Light curves are computed by integrating over the 
individual images. Repeating the procedure over observed 
wavelengths gives a time-dependent spectrum. 

To calculate fluid properties at each point on a ray, 
the spacetime coordinates of the geodesic are transformed 
from Boyer-Lindquist to the modified Kerr- S child coordi- 
nates used in the simulation (|McKinnevl l2006bh . Since the 
accretion flow is dynamic, light travel time delays along the 
geodesic are taken into account. Data from the sixteen near- 
est zone centres (eight on the simulation grid over two time 
steps) were interpolated to each point on the geodesic. 



2.3 Radiative Modeling 

Computing emission and absorption coefficients requires 
converting simulation fluid variables (pressure/internal en- 
ergy, mass density, and magnetic field strength) into an 
electron distribution function in physical units. The black 
hole mass sets the length and time scales, while the mass 
of the initial torus provides an independent scale and fixes 



2 Units with G = c = 1 are used except where noted otherwise. 
In these units, 1 M for M87 is c± 4/xas, ~ 10 15 cm, and ~ 9 hours. 



3 iDexter. Agol fe Fr agile (2009) found a radiative efficiency from 
energy lost to numerical reconne ction of e ~ . 1 for Sgr A* in a 
non-conservative simulation from Fragil e~et al.l (|2007l V 
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Figure 1. Ratio of magnetic to rest energy density (first and third panels) and specific enthalpy measured at infinity (second and fourth 
panels) for single time steps of the MBD (left two panels) and MBQ (right two panels) simulations. The azimuthally averaged data are 
shown in colour, while contours of b 2 / pc 2 = 0.1, 8tt and —hut = 1.0, 1.1 for the full range of azimuths are overplotted. Also overplotted 
in blue are lines of constant polar angle, 6 = 18°, 30°. For M87, 1 M ~ 4/ias. The jet is clearly identified in MBD, but only roughly 
corresponds to lines of constant 0. There is no persistent jet in the MBQ simulation. 



the accretion rate. The scalings are such that n and b 2 are 
proportional to the accretion rate. There is no consensus 
for the electron distribution or geometry responsible for the 
millimetre emission in M87. The presence of an extended 
jet at 7 mm indicates that the jet is at least compara- 
ble in luminosity to the disc at that wavelength. Since our 
model consists of a GRMHD simulation where a jet is pro- 
duced from accretion onto a black hole, we also include a 
disc component. The emission mechanism is taken to be 
entirely synchrotron radiation for both components. The 
disc electrons are assumed to be thermal, obeying a rela- 
tivistic Maxwell (Maxwell-Juttner) distribution. As in pre- 
vious models for Sgr A*, we assume a two-temperature flow 
with a constant ion-el ectron temperature ratio which is left 
as a free parameter ([G oldstom Quataert & IgumenshcheyJ 
120051 : iMoscibrodzka et al.l 120091 : iDexter et al.ll2010h . With 
this parameter, the electron temperature is computed from 
the gas pressure in the simulation using the ideal gas law. 
We use the approximate angle- dependent, u npolarized emis- 
sivity from lLeung, Gammie fc Nobld (|201lh . Previous mod- 
els of M87 have also included non-thermal disc emission, 
which could be importa nt for explaining the radio spectrum 
(jBroderick fe Loebll2009h . 



2.4 Jet Emission 

We assume that the jet emission is from non- thermal elec- 
trons, whose particle distribution is a power law in electron 
energy (Lorentz factor) with a constant index p between 
low- and high-energy cutoffs, 7min,max. The unpolarized syn- 
chrotron emission and absorption coeffi cients, j v and cl v , 
for this distribution in cgs units are (Legg & Westfold 1968; 
lYuan. Quataert fe Naravanll2003l : iDexterlboilh . 



ne 2 (p - \)v c 
= 2^(7^-7^) W 

X [G(x m in) — G(X max I 3 

ne 2 (p- l)(p + 2) 7 x ~^~ 2 



x [Ga(x m in) - Ga(x max ) ■ 



(1) 



(2) 



with Vc 3 tmc 9 > Xmin - max v l (7min,max^c), v is the emit- 
ted frequency, n is the particle density, B is the magnetic 
field strength, is the angle between the photon wave- vector 
and the magnetic field, and the power law synchrotron inte- 
uls are, 
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is the unpolarized synchrotron function and K a (y) is a mod- 
ified Bessel function. The functions G(x) and Ga(x) are tab- 
ulated for the desired values of p and interpolated for the 
emissivity calculation. Evaluating the double integrals for 
constructing t he interpolatio n tables can be sped up signif- 
icantly using (Westfold 1959), 
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Figure 2. Sample disc/jet spectra varying each parameter independently: 7 m in (top left), p (top middle), 77 (top right), m (bottom left, 
in units of 10 -5 ), and Ti/T e (bottom right). The default parameter values are 7 m i n = 50, p = 3.25, 77 = 0.05, rh = 10 x 10 -5 and 
Ti/T e = 10. 



to convert them to sums of single integrals. This form for the 
absorption coefficient is original, derived in the same fash- 
ion as for the emissivity in previous work. Both coefficients 
reduce to the standard approximate forms when x m in = 

and X m ax -> 00. 

This form of the emissivity is more accurate for M87 
than approximate forms commonly found in the literature 
that assume a frequency far from those corresponding to the 
cutoff Lorentz factors. In M87, the frequency corresponding 
to the low-energy Lorentz factor cutoff, for 7 m i n = 10 — 
100, is vo ~ 10 10_11 Hz, close to the frequencies of interest 
for mm-VLBI (230 GHz and 345 GHz). Taking the low- 
frequency cutoff into account broadens the spectrum and 
smoothes the turnover from optically thick to thin. 

The magnetic field strength everywhere is taken di- 
rectly from the simulation. For the jet emission, we need 
to calculate a non-thermal particle density. In magnetically- 
dominated regions such as the jet, the particle density and 
internal energy from the simulations are highly inaccurate 
due to the artificially enforced floor values used for numerical 
stability. Instead of using these compromised values, we scale 
the internal energy to the mag netic energy with a con stant 
of proportionality, 77 (cf. iBroderick &; McKinnevll201ol l: 



n n th rj 



8tt p - 



— (m e C 7min) 



(8) 



where m e is the electron mass and which implicitly assumes 
that all of the internal energy is in electrons, or equivalent ly 
that the thermal energy in all particles is negligible and that 
the positive charge carriers are thermal. 

For self-consistency, the rest energy of these non- 
thermal particles should still be less than the magnetic en- 
ergy density. This leads to the condition, 



m e p — 1 



87T7 m in 



(9) 



~ m e , P p — 2 

where m e , P correspond to leptons (baryons) producing the 
jet. The strictest condition on 77 is found by assuming a 
baryonic jet, in which case (for p = 3 — 3.5), 



< 25 7min 
' ~ 10 



(10) 



This inequality is satisfied in all our models as the maximum 
77 considered is 0.1. 



^nth 



(?) 



where u nt h is the non-thermal internal energy density and b 
is the magnetic field strength in cgs units. Then the particle 
density, n n th, is taken from, 



2.5 Disc/ Jet Boundary 

The GRMHD simulation consists of a solution for a single 
component fluid, neglecting electrons and dynamical effects 
on the local particle distribution. Defining a disc/jet solution 
then requires choosing a condition for the boundary between 
the two components. There are a few possibilities. 
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Figure 3. Viable spectra with Ti/T e = 10 (left), 5, (middle), 3 (right). The large red symbols are the observational data used as either 
upper limits (squares), data points (circles) or waveband-averaged data points (triangles). The small blue circles are the flux in the 
unresolved radio core of VLBI observations. 



The general structure of numerical simulations of black 
hole accretion flows consists of a dense, thick disc with scale 
height H/R ~ 0.2 — 0.3 centered on the equatorial plane, 
surrounded by a tenuous wind and then a polar jet. One 
straightforward method is to define the disc/jet boundary 
at a particular polar angle (~ 30°). 

The outflows in GRMHD simulations ar ise in 

magnetically-dom inated regions (|De Villiers et al.l [2005; 
McKinnev 2006bh. The second method is then to define the 
jet as anywhere that b 2 / pc 2 > /, where / ~ 0.1 — Sty de- 
pending how much of the "disc wind" is to be included in 
the jet. 

Finally, the disc/jet components can be separated as 
the boundary between inflow/outflow. In this case, unbound 
fluid elements constitute the jet {—hut > 1, where h is the 
specific enthalpy and ut is the time-component of the covari- 
ant fluid four- velocity, iDe Villiers et al.l [2005). Any combi- 
nation of these criteria can also be used. A comparison of 
the three criteria for an azimuthally-averaged time step of 
the simulation is shown in Figure [1] and as expected they all 
lead to roughly the same definition of the jet for MBD. The 
inflow/outflow and magnetically-dominated criteria are in 
particularly good agreement, while the constant slice in po- 
lar angle tends to cut out a significant amount of the jet base 
crucial for mm- VLBI images. We a lso show the same crite- 
ria ap plied to the simulation from iMcKinnev fc Blandfordl 
(2009) initialised with a l arge-scale quadrupolar magnetic 
field (MBQ in iDexter etHlfeoid ). where no persistent jet 
forms. Little of the polar region in this simulation would 
qualify as a jet by our criteria that fluid be magnetically- 
dominated or unbound. 

In the fiducial models discussed below, we use the con- 
dition b 2 1 pc 2 > 1, but our results are insensitive to the exact 
criterion used to separate the disc and jet components. 



Table 1. Fiducial Model Parameters 



Model 


T % /T e 


m(10- 4 ) 


Tmin 


P 


V 


DJ1 


3 


1 


50 


3.25 


0.05 


J2 


10 


1 


50 


3.50 


0.10 



black hole spin is also fixed at a/M = 0.92 since we only 
consider models from the MBD simulation. 

Each radiative model of M87 then has the follow- 
ing parameters: the ion-electron temperature ratio Ti/T e 
for the disc component, the accretion rate M and the 
disc/jet boundary selection criterion for both components, 
and 7min,max, P, and 7 7 for the jet component. We fix 7 max = 
10 5 throughout (see iBroderick &; Loebl 12009). leaving five 
free parameters. This is many more free parameters than 
in prior radi ative models of Sgr A* from GRMHD sim- 
ulations (iMoscibrodzka et al. 120091; IDexter, Agol &; Fragile 



2009; IDexter et alJl20ld : IShcherbakov. Penna fe McKinnev 
2010! ). Sample spectra showing the effects of independently 
varying the parameters are shown in Figure [21 The ion- 
electron temperature ratio Ti/T e (bottom right panel) de- 
termines the relative disc contribution. At Ti/T e — 10, 
the disc emission is negligible at all wavelengths while at 
Tj/ T e — 3, the best-fit ting ratio for this simulation for Sgr 
A* (|Dexter et alJl20ld ). the jet portion dominates except in 
the millimetre, where thermal disc emission produces a sub- 
millimetre bump. The normalization and peak frequency are 
affected by both j m in (for the jet portion, upper left panel) 
and m (bottom left panel). The fraction of magnetic en- 
ergy converted into non-thermal jet particles, 77, changes 
the normalization of the jet spectrum (top right panel). Fi- 
nally, p fixes the spectral slope between the millimetre and 
IR/optical emission. 



2.6 Model Parameters 

These are the required elements to define a radiative model 
of M87. Unlike in Sgr A*, the inclination of M87 is fairly well 
constrained from observe d superluminal motion in the jet 
(jHeinz fc Begelmadfl997h . implying i < 30°. We fix i = 25° 
in this work, but discuss the effect of varying it in § 16.31 The 



3 FIDUCIAL MODELS 

To identify viable models, we compute spectra from a 
single time step of the MBD simulation data over a 
grid spanning reasonable values of the various parameters: 
Ti/Te = (3,5,10), m = (1,2.2,4.6,10,22,46,100) x 10" 5 , 
p = (3,3.25,3.5,3.75), 7 min = (10,30,50,100,1000) and 
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r\ = (0.01, 0.02, 0.05, 0.10). The relative accretion rate is de- 
fined as rh = Mc 2 /L edd , where L ed d ~ 8 x 10 47 (M/6.4 x 
10 9 M©)ergs -1 is the Eddington luminosity. Although we 
only use a snapshot of the simulation, the shape of the spec- 
tra as well as the morphology of the millimetre wavelength 
images (see SJ) do not change significantly over the simula- 
tion, even as the total flux varies by a factor of a few (see 
33}. In £14.11 we discuss the effect of time variability on the 
predicted emission region sizes for mm-VLBI. 

The resulting spectra are fit to multi-wavelength 
obse rvations. We fit to average va lues in the opti- 
cal (ISparks, Biretta &; Macchettol Il996f ) and near infrared 
(jPerlman et al.ll200lL 12007) and to the measured values at 
3.3 mm and 1.3 mm (|Tan et al.ll2008h . The far infrared mea- 
surements are treated as upper limits due t o their possible 
conta mination by dust in the host galaxy (|Perlman et al.l 
l2007h . The radio data are also treated as upper limits, since 
as in Sgr A* the spatial extent of the simulation is too 
limited to model large-scale emission. The uncertainties are 
taken as 30% in all cases irrespective of measurement errors, 
since we are interested in finding qualitatively reasonable 
spectra rather than quantitatively constraining parameters. 

Examples of the many feasible models from the grid of 
spectra are shown in Figure [3] The lines are all spectra for 
which x 2 < 0.5, split up by ion-electron temperature ratio. 
The x 2 values are low because of our artificial inflation of 
the error bars, and the number of observational constraints 
is equal to the number of free parameters in our model. For 
both of these reasons we make no attempt to quantify a 
goodness of fit for the model spectra. 

These disc/jet models ar e much different from previous 
semi- analytic models of M87 ([Yuan 2000; Broderick &; Loebl 
2009). The jet and disc emission both peak in the millimetre, 
unlike truncated disc models where the disc emission peaks 
in the radio. The peak frequency of the disc spectrum and 
its normalization are fixed by Ti/T e and rh. Both parame- 
ters scale the normalization of the spectrum with its peak 
frequency, so that it is not possible to produce the observed 
radio emission from the GRMHD accretion flow. It could, 
however, be produced by thermal or non-thermal disc elec- 
trons at large radius outside of the simulation volume. Then 
the truncated disc would be most closely related to our mod- 
els with Ti/T e = 10, where the disc component contributes 
negligibly. The disc components of the model spectra with 
Ti/T e = 3 are similar to the time- averaged spectra from ax- 
isymmetric GRMHD s imulations ([Moscibrodzka et al.ll201ll : 
Hilb urn fe Lian3l201lh , except that we neglect the effects of 
Compton scattering. 

The jet spectrum has the additional degrees of freedom 
rj and 7 m in, which allow the possibility of a spectral peak 
in the radio. However, this is seemingly in conflict with 7 
mm VLB ! observations, which find ext ended mas scale jet 
emission (ljunor. Biretta &; Livid [l999h . In our jet models 
and those from Broderick fc Loebl (2009), the emission is ex- 
tended when optically thick to synchrotron self-absorption, 
and the existing VLBI observations suggest that the jet spec- 
trum is still rising in the millimetre. Moving the spectral 
peak of the jet component also requires extreme values of the 
relevant parameters, with ra > 10 -3 , 77 ~ 1 and 7min 00 1. 
For both reasons, we do not pursue these models further. 
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Figure 4. Total (solid), disc (dotted) and jet (dashed) spectra 
from the two fiducial models along with the observational data 
points. The error bars shown are from the measurements, whereas 
30% uncertainties were used for the fitting. In the model with 
Ti/T e =3 (left panel), the thermal disc electrons produce the mil- 
limetre emission. For small electron temperatures (right panel), 
the entire spectrum is from the non-thermal jet component. 



3.1 Viable Parameter Ranges 

We select the best-fitting spectra to the chosen data with 
Ti = 3, 10 as fiducial models representative of the set of 
viable possibilities. This parameter fixes the relative con- 
tribution of the disc. The parameters of the two fiducial 
models (DJ1 and J2) are listed in Table [1] and their spectra 
are plotted in Figure |4] showing the total spectra as well 
as their separate jet and disc components. The solid points 
are the total observed flux, while the open circles show 
only the fl ux from the unresolved radio core of VLBI ob- 
servations (|Paulinv-Toth et al.|[l98ll : ISpencer fe Junorlll986l : 
iBaath et al.lll992h . This is a better comparison for the ex- 
tremely small scale emission seen in our models, and gives 
reasonable agreement for frequencies > 10 10 Hz. 

Typical parameter values producing the millimetre 
emission in the fiducial models are n < 10 7 cm -3 , b ~ 10 
G, and T e c± 2 x 10 10 K. It is important to note that the 
observed emission depends on the electron temperature it- 
self, despite the common parametrization in terms of Ti/T e . 
Simulations with hotter ions would favor similar values of 
T e , but much larger values of Ti/T e . Typical non-thermal 
jet particle densities are n n th < 10 4 cm -3 . 



3.2 Jet Only or Disc/Jet 

With Ti/Te = 10, the disc emission is negligible at all wave- 
lengths. The jet spectrum peaks in the millimetre, with a 
power law tail extending to the IR/optical. Models with 
Ti/T e = 3 are still jet-dominated at low and high frequen- 
cies, but the thermal emission leads to a sub-millimetre 
bump. This not only leads to significant disc emission at 
frequencies of interest for mm-VLBI, but the thermal ab- 
sorption can also attenuate the jet emission, which for the 
inclination of M87 is mostly from the counter-jet (see SJ). 
The spectrum alone cannot distinguish between the jet or 
disc electrons producing the millimetre emission. 



4 IMAGE MORPHOLOGY 

Images of the fiducial models as well as their jet and disc 
components from the time step of the MBD simulation used 
for spectral fitting are shown in Figure [5] at 1.3 mm and 
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Figure 5. Images from total (top row), disc (middle row) and jet (bottom row) components for the two fiducial models (left two columns 
and right two columns) at 1.3 mm (first and third columns) and 0.87 mm (second and fourth columns). The colours are scaled linearly 
from blue to red to yellow to white, with a dynamic range of 60. The panel size is 100 x 100/ias. The images are taken from the same 
time step used for the spectra in Figure [4] and have been rotated 90° to roughly align with the position angle of the large-scale jet at 7 
mm. Both fiducial model images are crescents from the combined effects of light bending and Doppler beaming. The disc component of 
DJ1 is similar to previous models of Sgr A*. The jet component in J2 is dominated by the counter-jet, while in DJ1 the jet component 
comes from the forward jet since the counter-jet emission is absorbed by the disc. 



0.87 mm, the two wavelengths of interest for mm-VLBL The 
orientation angle of the M87 black hole spin axis projected 
on the plane of the sky can be reasonably assumed to align 
with the orientation of the 7 mm jet structure. This assumes 
that the jet is launched along the spin axis, and that the jet 
remains coherent on parsec scales. The images in Figure [5] 
have all been rotated to this favored orientation. 

Images of M87 dominated by thermal particles in the 
accretion flow ( DJ1) are nearly ide ntical to those of Sgr A* 
(cf. Figure 11 of lDexter et"aHl201Ch . Doppler beaming from 
the Keplerian velocity profile is significant even at the low 
expected inclination of M87 (i — 25° here), but weaker than 
for preferred inclinations of Sgr A* 60°). The image is a 
crescent from the combined effects of beaming, light bending 
and gravitational lensing. The emission region is in the inner 
porti on of the accretion flow ( r ~ 5M) near the mid-plane 
( see iMoscibrodzka et al.l 120091 : iDexter et all l2Q10h . 

The jet emission arises from near the pole at very small 
radius (r ~ 2 — 4M) , which falls mostly in the "shadow" re- 
gion of the image, connected with an observer at infinity by 
photons that intersect the black hole. The portions outside 
the shadow for lines of sight at the low inclinations consid- 
ered here (i < 40°) are in the counter- jet rather than the 
forward jet. Inside the shadow, a small fraction of the emit- 
ted photons reach a distant observer, strongly suppressing 
the brightness of the forward jet. When the disc emission is 



negligible (J2), the lines of sight to the counter-jet are opti- 
cally thin and the counter-jet dominates the jet image. The 
vertical velocity components are small at these small radii, 
so that the Doppler beaming is dominated by the azimuthal 
velocity component. This leads to the asymmetry in the jet 
component of the image. The jet image is still a crescent, 
although its intensity drops off sharply due to the strong 
radial dependence of the jet emissivity (j u ~ n n th^> 2 ~ r -4 
with B r>j r _1 near the spectral peak). When the disc emis- 
sion is significant, lines of sight to the counter-jet become 
optically thick, and the jet emission is dominated by the 
forward jet contribution (lower left two panels of Figure [5}. 

These are the first images of a jet launching region from 
a simulation. Images of the same time step are plotted in 
Figure \6\ for inclinations of 90° (edge-on, top left panel) 
to 0° (face-on, bottom right panel) in even steps of cosz, 
where i is the inclination angle. For edge-on viewing, there 
is no distinction between jet and counter-jet and the image 
is split into two bright lobes close to the black hole. The 
image is highly asymmetric due to strong Doppler beaming 
from helical motion in the jet base. Moving to lower inclina- 
tion, more and more of the forward jet falls into the shadow 
region where the path lengths are shorter and less of the 
orbital velocity is aligned with the line of sight. This causes 
the brightening of the counter-jet relative to the forward jet. 
At small inclinations, the forward jet is barely visible and 
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the image is dominated by the counter-jet. This is the op- 
posite of the behavior seen in extended jet emission where 
the material is ultra-relativistic and the radial velocity dom- 
inates. In that case, the forward- (counter-) jet is strongly 
Doppler boosted towards (away from) the observer at small 
inclinations. The counter-jet in M87 isn't visible even on 
mas scales at 7 m m. These jet image s are a lso much differ- 
ent from those in iBroderick fc Loebl (|2009j ). where the jet 
is launched farther from the black hole (r ~ 10M) than is 
found from the simulation assuming that n n th oc b 2 . When 
the jet is launched farther from the black hole, the forward 
jet dominates the images at 1.3 mm. For this reason, the jet 
images discussed here are sensitive to the uncertain mass 
loading in ultra-relativistic jets. Jet particle density profiles 
for M87 found from pair producti on calculations are qualita - 
tively similar to those used here (jMoscibrodzka et al.ll201ll ) . 

In both the disc and jet models the circular photon 
orbit produces a black hole shadow at 1.3 mm. At 0.87 mm, 
both images are more compact. The jet component in DJ1 is 
less significant, and the J 2 image is essentially just a single 
crescent from Doppler beaming of the jet material with a 
velocity component along the line of sight. 



4.1 Predictions for mm-VLBI 

The fir st mm-VLBI obse rvations of M87 were conducted re- 
cently (jFish et alJ l gOllh using telescopes in Arizona, Cali- 
fornia, and Hawaiijjbut the results are not yet available. We 
therefore make predictions for both current and future tele- 
scopes, assuming a geometry (inclination and orientation) 
based on larger scale jet observations. With sufficient base- 
line coverage, sensitivity and observing time it may be pos- 
sible to construct observational images in the future. How- 
ever, for now the small number of available telescopes re- 
quires that models be fit to observations in the u-v plane. 
We Fourier transform the images to visibilities, whose am- 
plitudes are shown along with the images for both mm-VLBI 
wavelengths in Figure [7] Possible locations of current mea- 
surements are shown as the green lines, and the visibility 
amplitudes are interpolated to those locations and plotted 
against baseline length in the top panel of Figure [8] In both 
fiducial models the current telescopes are at an orientation 
where the visibility amplitude decreases monotonically with 
baseline length. Fit with a symmetric Gaussian model, the 
source FWHM size is 43 (36) /xas for this single time step of 
model DJ1 (J2). 

However, the models are time variable and sensitive to 
the sky orientation. To make a more robust prediction, these 
effects are taken into account by fitting symmetric Gaussian 
models to visibilities from 120 images spaced evenly over the 
last ~ 2000 M of simulation time (~ 2 years for M87) with 
sky orientations of —75° ±30°. In both models, the inferred 
source size increases with the total 1.3 mm flux due to the 
emission region becoming optically thick. The inferred sizes 
(DJIfwhm and J2fwhm) can be expressed as power laws in 
flux as, 



4 Submillimeter Telescope Observatory, SMTO; Combined Array 
for Research in Millimeter Astronomy, CARMA; and James Clerk 
Maxwell Telescope, JCMT. 



DJIfwhm ~ (33 ±6) ^ — ^ M as 5 (H) 

J2fwhm ~ (39 ±8) ( ^j ^ ) ^ as ' ( 12 ) 

where the uncertainties are the la scatter both from the time 
variability and sk y orientation. For the observed total flux 
of 1.77Jy found by I Tan et al.l ([2008), the predicted sizes are 
both 44/xas. VLBI observations of Sgr A* have found smaller 
fluxes on micro- arcsecond scales than those from single dish 
observations, in which case the values of 33/xas and 39/xas 
for Uy may be more appropriate for M87. 

We can also interpolate the visibilities from the fiducial 
models to the baselines probed by telescopes to be added to 
the mm-VLBI array in the near future in Chik0 and Mex- 
ico0 and the results are shown in the bottom panel of Figure 
[8] In both cases, the black hole shadow is accessible to obser- 
vations between Mexico and Hawaii Q In the jet-dominated 
model, the shadow also appears on baselines between Chile 
and Mexico. Finally, the VLBI closure phase provides an ad- 
ditional constraint for any triangle of baselines, independent 
of that given by the visibility amplitude. Predicted closure 
phases for our two fiducial models as functions of position 
angle and the projection of the Hawaii/California/ Arizona 
baseline triangle on M87 are shown in Figure [9] The pre- 
dicted closure phase can deviate significantly from zero, 
indicative of asymmetric structure, especially in the jet- 
dominated model (J2). These results hold over the range 
of sky orientations and simulation time considered. 

The maximum brightness temperatures, kTb — 
I„c 2 /2v 2 in cgs units, where I v is the observed specific in- 
tensity, from the models are 3 x 10 10 K (DJ1) and 7 x 10 10 
K (J2). The brightness temperature tracks the image inten- 
sity, and its maximum is larger for a more compact image. 
In the disc model, the maximum brightness temperature is 
in good agreement with the typical temperature of the elec- 
trons producing the millimetre emission. 



5 VARIABILITY 

As discussed above, these radiative disc/jet models of M87 
are also time-dependent. Light curves from the fiducial mod- 
els are shown in Figure [TUl at 1.3 mm and 0.87 mm. The 
variability in the jet model is from one extended event, 
while the disc model light curve looks similar to that for the 
same simulation in the Sgr A* mod eling (lower right panel 
of Figure 6 from iDexter et al. I l2010h . This disc variability is 
caused by fluctuations in the particle density and magnetic 
field strength driven by MRI turbulence, and are strongly 
correlated with the accretion rate. Both models are consis- 
tent with the findi ng that M87 varies at about 1 Jyyr -1 
flSteppe et al. 1988). However, longer light curves would be 
required to make a statistically significant statement about 

5 Atacama Pathfinder Experiment, APEX; Atacama Submillime- 
ter Telescope Experiment, ASTE; Atacama Large Millimeter Ar- 
ray, ALMA. 

6 Large Millimeter Telescope, LMT. 

7 JCMT or Submillimeter Array, SMA. 
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Figure 6. Jet images at observer inclinations ranging from edge-on (top left panel) to face-on (bottom right panel). The colours are 
scaled linearly from blue to red to yellow to white, with a dynamic range of 60. The panel size is 100 x 100/ms. The arrow shows 
the orientation of the black hole spin axis, and the grid is the xy-plane. For edge-on viewing there is strong asymmetry from Doppler 
boosting and no distinction between the jet /counter-jet, but as the inclination decreases the image becomes circularly symmetric and 
the counter-jet becomes more prominent as more of the forward jet emission is captured by the black hole. 



the variability. This is particularly true in the jet model, 
where the variability is the result of one event. It is unclear 
whether this is due to a transient effect in the jet forma- 
tion/propagation or recurrent variable activity. 



6 DISCUSSION 

We have created the first radiative disc/jet models of M87 
based on GRMHD simulations. The jet is formed, colli- 
mated and accelerated self-consistently from an accretion 
flow within the ideal MHD approximation and the limita- 
tions of the numerical scheme. The disc/jet boundary is 
taken as a contour in the ratio of magnetic to rest mass 
energy density or in the particle specific enthalpy as mea- 
sured at infinity, with similar results in either case. The disc 
portion is modeled as thermal electrons with a constant ion- 
electron temperature ratio as in previous models of Sgr A*, 
while the internal energy of the power law electrons in the 
jet region is scaled as a fixed fraction of the magnetic energy 
density. The emission mechanism is synchrotron radiation in 
both cases. Two separate types of models can describe the 
high resolution spectrum of M87. In one class of models, 
the disc emission is negligible at all frequencies and the jet 
produces the entire spectrum. In the other, the jet produces 
the low and high frequency spectrum, while the millimetre 
emission is produced by thermal disc electrons. 

6.1 Favored Parameter Values 

The favored jet parameter combinations are 7mm = 30-100, 
p = 3.25 - 3.5 and r\ = 0.02 - 0.10. The favored av- 
erage accretion rate is rh = (1 — 2) x 10 -4 , or M = 
(1 — 2) x 10 _3 M©yr _1 . For the jet portion, the parameters 
m, 7] and 7min are partially degenerate. Ignoring important 
effects of optical depth and the cutoff frequency on the spec- 
trum, i^ max ~ ^M 2 7 min , while z/ max oc Vrh^ in . Jet spectra 
with 7min = 1000 peak in the near-IR or optical, inconsis- 
tent with VLBI observations of extended jet emission in the 



radio. Smaller rh is required at larger 7 m in to keep the spec- 
tral peak in the millimetre. The fraction of magnetic energy 
density in non-thermal electrons, 77, can be used to scale to 
the correct normalization for any 7 m in and rh as long as it 
remains less than unity and small enough so that the jet re- 
mains magnetically dominated (Equation [TO]). The spectral 
slope, p, is fixed by the near IR and optical observations, 
while the ion-electron temperature ratio Ti/T e determines 
whether the disc or jet produces the emission at millime- 
tre wavelengths. The favore d value for this sim ulation from 
Sgr A* models, Ti/T e = 3 (jPexter et al.ll2010h . leads to a 
disc-dominated millimetre image. 

6.2 Radiative Efficiency 

By ignoring radiation in the simulations, we have assumed 
that the accretion flow is radiatively inefficient. We can 
check this assumption by calculating the radiative efficiency, 
e = L/Mc 2 , where L ~ 2 x 10 42 ergs _1 is the bolomet- 
ric luminosity of the models. This gives e ~ (1 — 4) x 
10 -2 , so that the radiative models are radiatively ineffi- 
cient and marginally self- consistent. The effects of radia- 
tive cooling on the dynamical solution should be considered 
in future simulations. Using axisym metric GRMHD sim- 
ulations, Mos cibrodzka et al.l ([201 lh found lower accretion 
rates, M c± lO -4 M0yr -1 , and larger bolometric luminosi- 
ties, L rsj 10 44 ergs _1 , leading to unphysically large radia- 
tive efficiencies (~ 16, see their Table 1). Their models used 
a smaller black hole mass and T% = T e , which together could 
explain the small accretion rates found. The bolometric syn- 
chrotron luminosity L ~ nB 2 T?M 3 oc M 2 T 2 /M 3 . Using 
the factor ~ 2.1 smaller black hole mass and factor of 2 
larger electron temperature then requires an order of mag- 
nitude decrease in accretion rate to produce the observed 
flux, which explains the difference. Their models also include 
Compton scattering, which is neglected here. Extending our 
non-thermal jet spectrum to the X-ray does not significantly 
change its bolometric luminosity, but including Compton 
scattering would both increase the bolometric luminosity 
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Figure 7. Images (first and third columns) and corresponding visibility amplitudes (second and fourth columns) from the two fiducial 
models (rows) at 1.3 mm (first two columns) and 0.87 mm (right two columns). The colours are scaled linearly from blue to red to 
yellow to white, with a dynamic range of 60. The panel size is 100 x 100/ms. The baseline orientations available to current (green) and 
near future (white and pink) telescopes are overplotted as lines. At this orientation, set by the direction of the 7 mm jet, the black hole 
shadow is accessible to future observations on a Hawaii-Mexico baseline. 



and allow us to include a constrain t fro m the observed X-ray 
lumin osity (jWilson Yand 120021 ) . In iMoscibrodzka et all 
the spectral energy distribution from the simulation 
without radiative cooling peaks at v > 10 20 Hz. It is unclear 
whether Compton scattering would have such an important 
effect if included in our 3D models. 

We can crudely estimate the relative significance of 
Compton scattering for both models. In the jet model, the 
ratio of peak vF v for synchrotron and self- Compton compo- 
nents in a homogeneous blob with a power law distribution 
of el ectrons, ignoring a ll relativistic and viewing effects is 
(e.g., ISari fe Esinll200lh : 



____ aTni?7min7 



max 5 



(13) 



vF v (IC) 
i/F^Syn) 

where <jt is the Thompson cross-section and R is the size of 
the emitting blob. Scaling to reasonable parameters for our 
jet model gives, 



^/max 

10 5 



-1/2 



^(IC) / _R_\ ( U N /TminW 2 / 

vF u (Syn) ~ ' \5M J V10 7 cm" 3 J V 50 / V 

(14) 

This estimate is not self- consistent. We have scaled the elec- 
tron number density to the typical value for thermal disc 
electrons which provide the scattering optical depth, but 
then assumed they have a power law energy distribution. 
This is a conservative approach, since thermal electrons 
would have a much sharper high energy cutoff. 

For the disc model, we u se th e approximate local pre- 
scription given in lEsin et all (|l996h to calculate the Comp- 
tonization enhancement to the bolometric thermal syn- 
chrotron luminosity. This enhancement depends sensitively 
on the scattering optical depth and weakly on the seed 
photon energy. The scattering optical depth is taken to be 
r = ncrH, where a is the Thomson scattering cross section 
and H = (H/R)R is the disc height, where H/R ~ 0.2 



for the simulation considered here. Comptonization is esti- 
mated from the peak frequency of the synchrotron spectrum, 
Up ~ 10 1 1-12 Hz for the disc models. The estimate is mostly 
independent of v v within this range. Using typical values 
of n = 10 7 cm" 3 , R = 1 M, and T = 2 x 10 10 K, we find 
the relative Compton luminosity Lic/£s yn — 1. Numerically 
integrating the Compton cooling rate over the simulation 
domain gives a similar answer. However, this enhancement 
factor is extremely sensitive to the assumed scattering op- 
tical depth. For example, replacing the disc height H with 
radius R in the expression above gives Lc/Ls ~ 10, or a 
radiative efficiency e < 1, invalidating our assumption that 
radiation can be added after the fact without changing the 
(thermo) dynamics of the accretion flow. 

From both estimates, it is clear that Compton scatter- 
ing is an important mechanism and will significantly affect 
the observed X-ray emission. However, we estimate that it 
may only lead to an order unity correction in the bolomet- 
ric luminosity of both of our models. Including optically thin 
synchrotron and Compton cooling self-consistently is an im- 
portant goal for future simulations of M87, and may signifi- 
cantly affect the resulting radiative model. However, this is 
complicated by the uncertain ion-electron coupling. Includ- 
ing optically thin radiative cooling in standard single fluid 
simulations assumes perfect coupling between the ions and 
electrons, since both species are cooled equally. 

For completeness, we can also estimate the bolometric 
luminosity from bremsstrahlung: 



10 4 



( n V 

V10 7 cm" 3 / 



10 10 K 



1/2 



R 
10M 



erg s 



(15) 

a few orders of magnitude less than the synchrotron lumi- 
nosity. Including the large-scale accretion flow, with much 
smaller particle densities and temperatures but a much 
larger volume, will increase this estimate and both Compton 
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Figure 8. Model visibility amplitude vs. baseline length for cur- 
rent (top) and near future (bottom) telescopes. The visibility 
amplitude falls off monotonically with baseline length for cur- 
rent measurements and can be fit reasonably well with a circular 
Gaussian model. The black hole shadow appears as a local mini- 
mum in the visibility profile, and is accessible to a future baseline 
between Hawaii and Mexico. 

scattering and bremsstrahlung could be important emission 
mechanisms at X-ray energies. 

6.3 Images 

Assuming that the M87 jet propagates along the black 
hole spin axis, its orientation angle is constrained to be 
roughly — 75 =b 30° measured E of N. The inclination has 
b ee n estimated to be — 25° f rom the Lorentz factor of the 
jet (iHeinz Begelmanl Il997l : iBiretta. Sparks Macchettol 
Il999l ). For these parameters, we can make predictions for 
ongoing mm-VLBI observations. On current baselines, we 
predict that M87 should appear as a compact source sim- 
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Figure 9. Predicted closure phase from the two fiducial mod- 
els (left and right panels) at 1.3 mm as functions of the posi- 
tion angle of the black hole spin axis and the projection of the 
Hawaii/California/ Arizona baseline triangle on M87. For a large 
fraction of parameter space, the predicted closure phase differs 
significantly from zero, indicative of the asymmetric structure of 
the images. 



ilar to Sgr A*. The fiducial models have FWHM sizes of 
22 — 68/ias when fit with symmetric Gaussian models over 
the entire range of sky orientations and simulation time steps 
considered. For reasonable values of the total 1.3 mm flux, 
the sizes are largely in the range 33 — 44/ias, similar to the 
size of 37m£is found in the fi rst mm-VLBI observations of Sgr 
A* (jDoeleman et al.ll2008T ). 

The black hole shadow is accessible to future observa- 
tions on baselines between Mexico and Hawaii and possibly 
Mexico and Chile. Although the predicted Gaussian sizes for 
current telescopes are nearly identical for both fiducial mod- 
els, future observations with additional telescopes or epochs 
should be able to distinguish between the two. This is be- 
cause the 2D structure of the images is significantly different 
and future baselines will probe orientations where the pre- 
dictions differ substantially. Future mm-VLBI observations 
will also measure polarized emission on event horizon scales, 
and including polarization in the radiative transfer calcula- 
tions is a goal for future work. 

Both images are crescents, as in the case of Sgr A*. 
When the emission region is compact enough and the ve- 
locities are roughly Keplerian disc motion and helical jet 
motion, the resulting image is a crescent, qualitatively inde- 
pendent of the details of the physics in the innermost part 
of the accretion flow. The sky orientation (position angle) 
determines the orientations measured by VLBL Thus, our 
predictions for mm-VLBI with current and future telescopes 
are fairly robust as long as the assumed geometry is accurate 
and the millimetre emission region is compact (r < 10 M), 
despite the fact that the fiducial models chosen here are only 
representative of a range of possibilities for explaining the 
spectral properties of M87. The compactness requirement 
for the emission region is equivalent to assuming a steeply 
declining emissivity with radius, and a small optical depth 
(r < 1) to the emission region. 

At inclinations of 10° (40°), the images become less 
(more) Gaussian and more (less) circular. The crescent mor- 
phology and size predictions are less sensitive to changes in 
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Figure 10. Light curves from the two fiducial models (left and 
right panels) at 1.3 mm (top row) and 0.87 mm (bottom row). 
The disc variability is driven by MRI turbulence, as previously 
found in similar m odels of Sgr A* (jPexter. Agol fc Fr agile 2009; 
iDexter et aDl2010h . 

inclination within this range than to uncertainties in the 
model (DJ1 or 32) or position angle. 

We have assumed the new mass estimate of 6. 4 x 1 9 M 
(jGebhardt fe Thomas! [20091 : 1 Gebhardt et al.ll201lh through- 
out. Changing the black hole mass changes the favored pa- 
rameter values, but does not substantively change the types 
of feasible models or other major results. The exception is 
in the size predictions for mm-VLBI, which essentially scale 
with mass. If the black hole mass in M87 is 3 x 1O 9 M 
(jMarconi et al.lll997T ). then the predicted sizes would be a 
factor ~ 2 smaller. This means that mm-VLBI can provide 
a model-dependent test of the black hole mass - if a size 
smaller than ~ 22/xas is found, the smallest for our models 
over the range of simulation time and valid position angles, 
it may indicate a smaller black hole mass. Conversely, a rel- 
atively large measured size would favor the larger black hole 
mass estimate, although it could also mean that the emis- 
sion at 1.3 mm is more extended than found from current 
GRMHD simulations. 

6.4 Uncertainties 

There are many uncertainties in this analysis. The as- 
sumption that the internal energy in non-thermal particles 
scales with magnetic field energy densit y may be reason- 
able (e.g., iBroderick &; McKinnevI l2010h , but it is made 
out of necessity. Both the internal energy (pressure) and 
mass density from the simulation are set by the artificial 
numerical floor required for code stability when b 2 / pc 2 ^> 
1, the region of interest for jet launching. Radiation is 
added in post-processing, despite the fact that the models 
are found to be marginally radiatively efficient. Both syn- 
chrotron and Compton scattering have been found to be 
dynamically important in axisymmetric simulations of M87 
( Mos cibrodzka et ajj [201l| ) with similar parameters to our 



disc/jet model (DJ1). It will be possible to include those 
forms of cooling in future simulations, but a method for 
evolving the non-thermal particle density self-consistently is 
much more difficult. 

Previous semi-analytic disc/jet models of M87 have in- 
voked a truncated disc, with a low, constant particle den- 
sity and magnetic field strength throughout the inner disc. 
This configuration can produce the observed radio emis- 
sion. We find that the magnetic field strength in MBD and 
other GRMHD simulations falls off with radius (oc r -1-2 
for MBD). Producing the observed flux then requires that 
the disc emission peaks at millimetre wavelengths. Although 
none of our models can explain the observed radio emission, 
the agreement is much better if the models are fit to only 
the unresolved emission from the core of r adio images (cf. 
iRevnolds etaDll99d : iDi Matteo et al.ll2003h . 

Both the disc and jet models used here are unlikely to 
be valid outside the innermost radii. The accretion flow so- 
lution is probably reliable only out to r ~ 10 M, but the 
disc emission region is so compact that this does not af- 
fect the results. However, there are so few studies of jets 
launched from GRMHD simulations that it is unclear what 
the domain of validity is. This situation will improve in the 
near future, as numerical techniques for propagating large- 
scale jets improve and simulations are able to be evolved 
for longer physical times due to the increase in available 
computational resources. The jet solution in a new sim- 
ulation similar to MBD but with much higher resolution 
(272x128x256), a larger radial boundary (26000 M) and 
a longer duration (20000 M) is nearly identical to that of 
MBD within r < 8 M (B oc r -1 ), but has a nearly flat 
radial magnetic field strength profile outside of that. Pre- 
liminary images based on a single time step from this sim- 
ulation display extended (mas) structure in the radio, but 
give nearly identical predictions for mm-VLBI of M87, in- 
dicating that the low resolution used in MBD is unlikely 
to invalidate the results of this paper. Models from this 
and other high resolution, geometrically thic k simulations 
(Mc Kinney, Tchekhovskov & Blandfordl f2012) will be con- 
sidered in a subsequent publication. 



7 SUMMARY 

We have constructed the first radiative disc/jet models of 
M87 based on a general relativistic MHD simulation of a 
black hole accretion flow. The main results of this study 
are: 

(i) The spectral energy distribution of the core of M87 
can be explained with jet or disc/jet models. 

(ii) In both types of models the images at 1.3 mm are 
crescents from the combination of gravitational lensing and 
Doppler beaming of the compact emission region. 

(iii) The jet emission is produced in the counter-jet near 
the pole in the immediate vicinity of the black hole (r < 4 
M), while the disc emission is produced near the midplane 
in the inner radii r c± 5 M as found in previous models of 
Sgr A*. 

(iv) For the favored viewing geometry of M87 based on 
analyses of the large-scale jet, we predict a Gaussian source 
structure with a size of 33 — 44/xas (~ 5 Schwarzschild radii) 
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for observations with current mm-VLBI telescopes. The in- 
ferred size should increase slowly with increasing flux. The 
black hole shadow, direct evidence for an event horizon, may 
be detected on future baselines between Hawaii ( JCMT) and 
Chile ( ALMA/ APEX/ ASTE). 

(v) The two types of models can be distinguished with 
mm-VLBI observations on baselines including telescopes in 
Mexico (LMT) and Chile ( ALMA/ APEX/ ASTE), or by 
measuring the change in size between epochs of different 
total 1.3 mm flux. 
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